function vmd_Uspectra(fbase, fnum)
%
%   vmd_Uspectra(fbase, fnum)
%
%   Reads file with density distribution.
%   Creates two text files with specral data.
%
%   Input:
%        fbase   base of file number (runname)
%        fnum    file number
%
%   Domain information and colormap limits are edited in the file.

%---------------------------------------------------------

% system size

%jx=48; jy=48; jz=16; jv=256; jdir=6;
jx=48; jy=1; jz=16; jv=256; jdir=2;

% physics paramters

Lx=450; Ly=450; kz=0.34; nz=1;


%---------------------------------------------------------

  dkz = kz/nz;
  dkx = 2*pi / Lx;
  dky = 2*pi / Ly;

  ixc = jx/2+1;
  iyc = min(jy/2+1, jy);

  fn = num2str(fnum,'%04d');


  %-- read velocity components

  fname = [fbase, '.ux.', fn];
  fid  = fopen(fname, 'rb');
  ux   = fread(fid, jx*jy*jz*jdir, 'double'); fclose(fid);
  ux   = reshape(ux, jz,jy,jx,jdir);


  fname = [fbase, '.uy.', fn];
  fid  = fopen(fname, 'rb');
  uy   = fread(fid, jx*jy*jz*jdir, 'double'); fclose(fid);
  uy   = reshape(uy, jz,jy,jx,jdir);


  fname = [fbase, '.uz.', fn];
  fid  = fopen(fname, 'rb');
  uz   = fread(fid, jx*jy*jz*jdir, 'double'); fclose(fid);
  uz   = reshape(uz, jz,jy,jx,jdir);



  %-- output files --

  fname = [fbase, '.UmatFTx.', fn, '.txt'];
  fidx = fopen(fname, 'wt');

  fprintf(fidx,'%% z-averaged x-spectra of U-field\n');
  fprintf(fidx,'%% 1.kx   2.<|Ux(kx)|^2>   3.<|Uy(kx)|^2>   4.<|Uz(kx)|^2>\n')

  fname = [fbase, '.UmatFTz.', fn, '.txt'];
  fidz = fopen(fname, 'wt');

  fprintf(fidz,'%% z-spectrum of Uz along center line\n');
  fprintf(fidz,'%% 1.kz   2.<|Uz(kz)|^2>\n')


  %-- loop over directions --

  for idir = 1:jdir

    fprintf(fidx, '\n\n%% Direction = %d\n\n', idir);
    fprintf(fidz, '\n\n%% Direction = %d\n\n', idir);

    Ux   = ux(:,:,:,idir);
    Uy   = uy(:,:,:,idir);
    Uz   = uz(:,:,:,idir);

    ux2=sum(sum(sum(Ux.*Ux)))/(jx*jy*jz);
    uy2=sum(sum(sum(Uy.*Uy)))/(jx*jy*jz);
    uz2=sum(sum(sum(Uz.*Uz)))/(jx*jy*jz);

    fprintf("Direction %1d mean squares (x,y,z): %14.6e %14.6e %14.6e\n",...
	    idir, ux2, uy2, uz2);    


    %-- x-spectra

    Uk   = fftn( Ux(:,iyc,:) )/(jx*jz);
    Uksq = Uk.*conj(Uk);
    Uxsq = sum(Uksq,1)/jz;

    Uk   = fftn( Uy(:,iyc,:) )/(jx*jz);
    Uksq = Uk.*conj(Uk);
    Uysq = sum(Uksq,1)/jz;

    Uk   = fftn( Uz(:,iyc,:) )/(jx*jz);
    Uksq = Uk.*conj(Uk);
    Uzsq = sum(Uksq,1)/jz;

    for i=1:jx/2
       k=(i-1)*dkx;
       fprintf(fidx, "%14.6e %14.6e %14.6e %14.6e\n", k, Uxsq(i), Uysq(i), Uzsq(i));
    end


    %-- z-spectra 

    U = Uz(:,iyc,ixc);
    U = reshape(U,jz,1);

    Uk = fft(U)/(jz);
    Uksq = Uk.*conj(Uk);

    for i=1:jz/2
      k=(i-1)*dkz;
      fprintf(fidz, "%14.6e %14.6e\n", k, Uksq(i));
    end

  end

  %-- end of loop over directions --
  
  fclose(fidx);
  fclose(fidz);

end
